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Abstract 

We present a Monte Carlo study of external surface anchoring in nematic cells with partially dis- 
ordered solid substrates, as well as of intrinsic anchoring at free nematic interfaces. The simulations 
are based on the simple hexagonal lattice model with a spatially anisotropic intermolecular po- 
tential. We estimate the corresponding extrapolation length b by imposing an elastic deformation 
in a hybrid cell-like nematic sample. Our estimates for b increase with increasing surface disorder 
and are essentially temperature-independent. Experimental values of b are approached only when 
both the coupling of nematic molecules with the substrate and the anisotropy of nematic-nematic 
interactions are weak. 

PACS numbers: 61.30. Cz, 61.30.Gd 
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I. INTRODUCTION 



In confined nematic liquid crystals with a large surface-to-volume ratio the aligning ef- 
fects of the confining surfaces are of great importance in determining the equilibrium director 
configuration [1]. There are two major contributions to these surface aligning effects, the first 
one originating from direct interactions between nematic molecules and the solid substrate 
(external anchoring), while the second one is due to incomplete anisotropic nematic-nematic 
interactions in the vicinity of the sample surface (intrinsic anchoring). An understanding 
of these confinement related aligning mechanisms is of great importance not only from the 
fundamental point of view, but also, e.g., for the design and construction of liquid crystal- 
based optical devices. Within phenomenological approaches anchoring effects are usually 
characterized by two parameters: the preferred molecular alignment direction at the sample 
surface (the easy axis) and the free energy coefficient W penalizing any deviation from this 
direction (the anchoring strength) [2]. Here the W coefficient depends on S, the nematic 
order parameter. This dependence seems to be strongly related to the specific properties 
of a given interface: for example, for a system of hard rods confined between hard walls 
one finds W oc S [3], while in experiments measuring anchoring inside polycarbonate mem- 
branes even W oc S* 4 can be obtained [4]. Given, moreover, K oc S 2 [5, 6] (K denoting 
the Frank elastic constant), with decreasing S the extrapolation length b = K/W [5] may 
either increase (as seen experimentally in thermotropics [4]), or decrease (as obtained from 
simulations with hard particles [7]). In the bulk, the value of S is primarily determined by 
temperature, while close to an interface it may also be affected by incomplete bulk interac- 
tions, as well as by (dis)ordering effects of the (possibly rough) confining substrate. In the 
latter case W is often assumed to be simply proportional to WS S(0), where W represents 
the surface coupling constant, So the surface-imposed value of S (determined by substrate 
roughness), and S'(O) its actual surface value [8]. 

While both the anchoring strength and the easy axis can be determined experimen- 
tally [9], it is also possible to deduce them from simulations based on pairwise intermolec- 
ular interactions. For example, intrinsic anchoring has been analyzed in a pseudomolecular 
continuum approach with ellipsoidal molecules [10], using a lattice approximation in the 
zero-temperature limit [11], or analyzing different types of interfaces in Gay-Berne sys- 
tems [12, 13, 14]. Furthermore, surface anchoring strength has also been measured in a 
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system of hard ellipsoids in contact with a hard wall [7, 15]. The anchoring strength re- 
ported in most analyses shows that anchoring is rather strong and that the corresponding 
extrapolation length is of the order of a few molecular dimensions, while its experimental 
values are typically above 100 nm [9]. On the other hand, a recent analysis of external 
anchoring — also based on a lattice model, but for nonzero temperature — was presented in 
Ref. [16], showing that the extrapolation length can increase significantly when the nematic- 
to-isotropic (NI) transition is approached, as also seen experimentally [4]. Moreover, a lattice 
gas approach has been adapted recently to study nematic interfaces [17]. 

Motivated by these developments, in this paper we extend the analyses performed in 
Refs. [11, 16] to nonzero temperatures (employing Monte Carlo simulations) and repeat 
the measurement of the extrapolation length both for external and intrinsic anchoring, the 
former in presence of rough solid substrates and the latter for a free nematic interface. The 
term "rough" here refers to flat but partially disordered substrates. We will first briefly recall 
the features of the lattice model used in Ref. [11], then discuss the modifications needed to 
perform the present analysis, and finally, present and discuss the results. 



II. SIMULATION MODEL AND ANCHORING MEASUREMENT 

We model a liquid crystalline slab of thickness d, using a modification of the well-known 
Lebwohl-Lasher (LL) model [18]. In the LL model elongated nematic molecules are rep- 
resented by freely-rotating spin-like particles that are attached to lattice points of a cubic 
lattice. Since the pairwise interaction energy used in the LL model is spatially isotropic 
and thus does not depend on the relative position of the particles, it cannot produce any 
orienting effects at a free nematic surface and is therefore not suitable for studies of intrinsic 
anchoring. Therefore, we choose a more general potential to model the interparticle interac- 
tion energy resulting from anisotropic van der Waals forces. These forces can give rise to the 
nematic phase [19] even without hard rod repulsion. For two neighboring nematic particles 
i and j with orientations given by unit vectors Uj and Uj, separated by r^, we model the 
interaction potential as [20]: 

1 2 

(1) 



Uij = -e 



Ui -uj -3z/(u,- ^)(u, 

I ij I ij j 
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with e > 0, which for v = reduces to the standard LL potential, while for v — 1 one has 
the induced dipole-induced dipole interaction. The spatial anisotropy parameter v enables 
us to continuously vary the relative importance of the spatially anisotropic (r^-dependent) 
contribution to the interaction law (1). Although the range of van der Waals potential 
is proportional to r^ 6 , for computational efficiency we consider only interactions between 
nearest neighbors. Thereby the intrinsic anchoring energy W and the elastic constant K are 
underestimated, but we expect the errors in the estimation of W and of the extrapolation 
length K/W not to exceed 20% and 30%, respectively. 

In a nematic slab, boundary conditions are typically fixed by the interaction with solid 
walls, or, alternatively, through orienting effects near free nematic interfaces. In the simu- 
lation, each of the solid walls is represented by a layer of fixed particles (]pj| = 1), either 
all perfectly aligned or somewhat disordered with some residual order. The nematic-wall 
interaction is modeled via 

^. = - es [p,-u/, (2) 

as promoted, e.g., by short-range steric forces, where, again, pj and u^- are nearest neighbors 
on the lattice. A dimensionless anchoring strength parameter can be defined as w = e s /e. 
In the case of a perfectly aligned surface one has p; = II, i.e., all particles are aligned along 
the easy axis II, while for disordered surfaces II represents the average orientation of p«. 
On the other hand, a free nematic interface is simply modeled through the missing-neighbor 
effect. 

Instead of the cubic lattice used in the LL model, in present simulations we use a simple 
hexagonal lattice to model the liquid crystal, which is necessary to avoid unphysical bulk easy 
axes [11] for v ^ 0. These axes arise as a direct consequence of using a lattice approximation 
for modeling an anisotropic liquid and are present for any spatially anisotropic potential in 
the cubic lattice. They are absent, however, in the hexagonal lattice, provided that the spins 
Uj are assumed to be two-dimensional vectors confined to hexagonal planes. The geometry 
of our sample is shown in Fig. 1. The z axis is normal to the confining interfaces and the 
hexagonal planes are parallel to the xz plane. We stress that our analysis must be restricted 
to low values of the anisotropy parameter v in order to avoid solid-like periodic director 
profile solutions that are stable for v > 0.3 [11]. 

The strength of any anchoring can be measured by imposing an elastic distortion so that 
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FIG. 1: Sample geometry, hexagonal lattice and the three sublattices (squares, circles, and trian- 
gles): (a) xz and (b) xy cross sections. The tilt angle <f> = 4>{z) is measured with respect to the 
z-axis, the sample normal. Boundary conditions at z = and z = d are homeotropic and planar, 
respectively. 



the average surface molecular orientation deviates from the easy axis denned by the anchor- 
ing. The magnitude of this deviation can then be used to estimate the anchoring strength 
and the corresponding extrapolation length [5]. The elastic distortion in a nematic slab can 
be imposed either by applying a magnetic field whose orientation must not coincide with the 
direction of the easy axis, or by antagonistic anchoring conditions at opposing surfaces. In 
the zero-temperature analysis of Ref. [11] the magnetic field method was used. On the other 
hand, for nonzero temperatures a strong enough magnetic field can enhance the degree of 
nematic order and even shift the NI phase transition, which can present additional difficul- 
ties in interpreting the results. Therefore, we decided not to use the magnetic field approach; 
instead, we consider a hybrid cell-like sample with antagonistic boundary conditions. 

Consider the v = case with pure external anchoring first. The left (z = 0) interface 
- where the anchoring will be measured — is chosen to be a solid wall, either perfectly 
ordered or somewhat disordered, promoting homeotropic alignment through weak external 
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anchoring (w = 0.1). At the same time, the right (z = d) interface is also taken to be a 
solid wall, yet with perfect planar alignment and strong external anchoring (w — 1). For 
the case v ^ 0, we replace the solid wall at z = with a free nematic interface, enabling us 
to study intrinsic anchoring alone. Note that for v < 0.3 the easy axis of the free interface 
is homeotropic [11], thus providing the same confinement type as with external anchoring. 
As a result, in both cases a combined bend and splay elastic deformation is expected to 
appear in the sample. The deformation should be present as long as the sample thickness 
d exceeds d c = \(K/W) — (K/W) d \, where (K/W) and {K/W) d are the extrapolation 
lengths corresponding to the effective anchoring on the left and the right wall, respectively, 
while for d < d c a uniform director structure is observed [21]. In a hybrid cell close to NI 
transition a non-bent biaxial structure is also possible, consisting of two strata of uniform 
alignment, separated by a biaxially ordered layer [22]. Recall that unlike in the original 
LL model in the present study spins Uj are two-dimensional vectors and hence unable to 
reproduce biaxiality. 

For hybrid boundary conditions in the one-constant approximation the Frank elastic 
theory (assuming a constant degree of nematic order throughout the slab) predicts a perfectly 
linear director tilt angle profile <f>(z), where is measured, e.g., with respect to the slab 
normal z. The torque balance condition at, e.g., the left surface (z = 0) can be written 
as (d(f)/dz) = \(W/K)q sin 20(0) and enables us to deduce (K/W) from the measured 
values of (d<fi/dz) and 0(0). Note that for small 0(0) (i.e., strong enough anchoring) the 
above condition simplifies to (d<j)/dz)o = (W/K) (f)(0) and allows us to determine (K/W) 
simply by extrapolating the profile <f>(z) graphically across the sample boundary to = 
corresponding to the homeotropic easy axis. Note also that if the degree of nematic order 
is subject to variations — which is usually the case near sample boundaries in a layer of 
thickness £ (£ denoting the nematic correlation length) — the profile (f>(z) may deviate from 
the above predicted linear behavior. In this case the extrapolation of the profile towards the 
surface must be performed from far enough in the bulk where the order parameter profile is 
constant, that is from £ < z < d — £. Note also that whenever (K/W) approaches d (for 
weak anchoring or in a thin sample), (d<p/dz)o and, consequently, (K/W)q are accompanied 
by a significant systematic error. 

Our Monte Carlo (MC) simulations are now performed as follows. The size of the simula- 
tion box size was set to 48 3 for a total of 105984 nematic particles, excluding the boundary 
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FIG. 2: Hybrid cell with smooth substrates and v = (external anchoring at 2 = 0): director (a) 
and order parameter (b) profiles for T* = 1.50, 1.42, 1.34, and 1.26 (diamonds, squares, triangles, 



and circles, respectively); the bulk NI transition temperature approaches Tfoj ~ 1.52. 



particles Pi in the layers at z = and z = d. In case of disordered substrates parti- 
cle orientations in the confining layers are generated following a probability distribution 
f(<f>) oc exp (— V cos 2 <p), which gives a homeotropic easy axis ((f) = 0) and a "P-dependent 
degree of order So. For example, V — > corresponds to a distribution close to isotropic, 
while V — > oo yields a perfectly aligned substrate. Then, the hexagonal lattice is divided 
into three sublattices, as shown in Fig. 1, ensuring that the bonds between neighboring 
particles on the lattice never connect two particles from the same sublattice. Considering 
the simple hexagonal lattice as tripartite enables us to vectorize the simulation algorithm, 
which provides a significant speed-up in calculations. Further, in the x and y directions peri- 
odic boundary conditions are assumed. We start either from a random configuration in two 
dimensions (recall that are restricted to lie in hexagonal planes), or from an equilibrated 
configuration at a temperature slightly higher than the simulated one, if this is available. 



Then we apply the standard Metropolis algorithm [23]. For our vectorized algorithm to 
work correctly, in each MC cycle we first attempt (and accept /reject) trial moves involving 
particles in the first sublattice and only then proceed to the second one, and after that 
to the third one. In generating a new trial configuration each time only a single particle 
is involved. We have typically performed 10 5 MC cycles for equilibration, followed by 10 5 
production cycles to accumulate averages of interest. 

We measure the extrapolation length by analyzing the director profile 4>(z). The <p(z) 
dependence is calculated by accumulating the two independent components of the two- 
dimensional ordering matrix Q a p{z) = 2( u ? u i)z — S a /3, where a and (3 can be either x or z, 
and uf represents the a component of the unit vector Uj. The average (...) z is performed 
both over all particles in the layer centered at z and over the production MC cycles. Then, 
the averaged ordering matrix is diagonalized and the positive eigenvalue is identified as the 
two-dimensional scalar order parameter s(z). Accordingly, the corresponding eigenvector is 
the director. Note that in Ref. [16] the extrapolation length in the LL model with three- 
dimensional spins was studied, using a different method to obtain the director profile. Rather 
than computing the ordering matrix, the polar angle averaged over all particles in a layer and 
over production MC cycles was computed; i.e., no information was kept about the azimuthal 
orientation of the particles. This method overestimates the actual average director tilt angle 
for small values of this angle. The strong temperature dependence of the extrapolation 
length found in Ref. [16] is most likely an artifact of the incorrect method used in that 
paper. 



III. EXTERNAL ANCHORING 



We analyze the spatially isotropic intermolecular interaction first by setting v = in 
Eq. (1). Such a model is now similar to the standard LL model, yet it is characterized by a 
different coordination number (eight in the former model, six in the latter) and the dimen- 
sionality of nematic spins (two in the former model, three in the latter). These distinctions 
lead to a different balance between the decrease of internal energy and loss of orientational 
entropy upon going from the isotropic into the nematic phase, and thus shift the NI transi- 
tion to a temperature higher than in the LL model. Monitoring temperature scans of internal 
energy, specific heat, and the order parameter s (lowercase s will be used throughout the 
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FIG. 3: Hybrid cell with a rough homeotropic substrate at z = 0: director (a) and order param- 
eter (b) profiles for v = 0, T* = 1.26, and different degrees of substrate roughness so- Circles, 
triangles, squares, and diamonds correspond to sq = 1, so ~ 0.78, sq « 0.41, and so « 0.15, 
respectively. 

text for the two-dimensional order parameter), for u = we estimate the dimensionless 
transition temperature in a 48 3 bulk sample (with full periodic boundary conditions) to be 
Tni = ksTNi/t ~ 1.52 ± 0.01. Here a dimensionless temperature scale, T* = ksT/e, has 
been introduced. 

Consider the director and order parameter profiles [(f>(z) and s(z), respectively] in a hybrid 
nematic cell where both substrates are smooth and hence impose perfect nematic order with 
So = 1. Recall that for v — one is dealing exclusively with external anchoring. Fig. 2(a) 
shows <p(z) profiles for different temperatures. One can readily observe that far enough 
from the NI transition the profiles approach a linear function (predicted also from Frank 
elastic theory [5]), with minor changes in slope \d<j)/dz\ only close to substrates where the 
degree of order can vary. In particular, wherever s(z) exceeds its bulk value s&, the nematic 
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becomes more difficult to deform, which results in a reduction in slope (and vice versa). 
Moreover, note that far enough from the NI transition [in our simulations for values of 
reduced temperature r = (T^ I — T*)/T^ I > 0.07] the 4>(z) profiles are essentially insensitive 
to changing temperature. Given w — 0.1, as chosen above for the homeotropically anchored 
smooth left wall, the extrapolation length is estimated as b = 11.6a(l ± 7%), where a is 
the lattice spacing. Significant deviations from the linear (f>(z) profile can be observed only 
close to T^j [for r < 0.02], when the nematic far enough from the walls melts and thereby 
avoids elastic distortion — see Fig. 2 (b). Then molecular alignment becomes homeotropic 
in the vicinity of the left surface, followed by a region of (nearly) isotropic liquid in the 
slab center, and by a region of planar alignment close to the right surface. Note that one 
can avoid this nematic "meltdown" by reducing the deformation strength, e.g., by setting 
< (f)(d) <C 7r/2 or by reducing w at the z = d wall. In this case b is observed to be 
essentially T*-independent even up to r ks 0.02. Note, however, that thereby the immediate 
vicinity of the NI transition actually has not been probed: for a realistic liquid crystal with 
T NI ~ 300 K, r 0.02 corresponds to temperatures as much as 6 K below the transition. 
Alternatively, the deformation strength can be reduced also by significantly increasing the 
system size. Further, in Fig. 2 (b) one can always observe an increase of s when approaching 
the substrate at z — d promoting planar alignment. Similarly, s decreases close to the z = 
substrate because of a weaker nematic-substrate coupling (w — 0.1 as opposed to w — 1 for 
the z — d substrate). 

Turn now to cases with a rough homeotropic substrate, simulated by generating an en- 
semble of fixed particles pi at z = with s < 1. Qualitatively, the behavior of <f>(z) and 
s(z) does not change, including the insensitivity of b to variations of T* in the temperature 
range presently accessible. However, as shown in Fig. 3 for T* = 1.26, a smaller sq reflects in 
a reduced anchoring strength at z = and hence in a reduced slope \d(f>/dz\, which, in turn, 
results in a larger b: for s ^ 0.78 we find b = 14.9a(l±7%), for s ^ 0.41, b = 31a(l±7%), 
and for s rts 0.15, b = 50a(l±48%) (for a summary see Table I). Note that in the latter case 
b already approaches the sample thickness d and is already close to violating the stability 
condition d ^ b for the deformed director structure [21]. Moreover, assuming K oc si and 
W oc So s (0), the extrapolation length should scale as b oc sl/s(0)s . The full quantitative 
agreement, however, turns out to be rather poor (even within error bars). 

Note that in the vicinity of fairly disordered substrates the nematic becomes elastically 
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TABLE I: Extrapolation lengths b measured in units of lattice spacing a: substrate roughness- 
dependence for external anchoring (y = 0) and values for intrinsic anchoring (y ^ 0), together 
with the corresponding dimensionless bulk NI phase transition temperatures Tfcj. The values are 
given for T* < 1.42 [y = 0), T* < 1.265 (y = 0.05), and T* < 1.13 {v = 0.1), i.e., for r > 0.07, 
where b is essentially temperature-independent. 



V 


1 NI 


anchoring 


«0 


b (a) 





1.52 ±0.01 


external 


1.0 


11.6(1 ±7%) 








0.78 


14.9(1 ±7%) 








0.41 


31 (1 ± 7%) 








0.15 


50 (1 ± 48%) 


0.05 


1.36 ±0.01 


intrinsic 




16.3(1 ± 13%) 


0.1 


1.21 ±0.01 


intrinsic 




4.6 (1 ± 10%) 



softer than in the bulk, which increases the local slope of the director profile (f>(z) in com- 
parison with the bulk. As the NI transition is approached, the thickness of this disordered 
layer (~ £) starts to grow, thereby indirectly affecting the bulk slope of <f>(z), which — in 
principle — could result in a more pronounced temperature dependence for b that is extrap- 
olated from the bulk. Note, however, that a significantly thicker sample than the present 
one is required to actually observe this scenario. 

We also studied the influence of the hexagonal lattice structure and the confinement of 
Uj to the xz plane on our results, by analyzing external anchoring in the original LL model 
where the spins Uj are three-dimensional vectors. Both the LL and hexagonal lattice models 
yield the same qualitative behavior: a nearly temperature-independent b which increases in 
value as the surface roughness is increased. 

IV. INTRINSIC ANCHORING 

Setting v ^ 0, intrinsic anchoring appears at the interfaces, in addition to the external 
contribution present already for v — 0. Note that changing v affects not only anchoring, 
but also changes the elastic softness of the nematic [24]: for v not too large [y < 0.78), the 
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FIG. 4: Nematic slab with a single free interface at z = and v = 0.05 (homeotropic intrinsic 
anchoring). Same profiles as in Fig. 2, yet for rescaled values of T*: T* = 1.335, 1.265, 1.195, and 
1.125 (diamonds, squares, triangles, and circles, respectively). T^j m 1.36. 

Frank elastic constant decreases with increasing is. Consequently, upon increasing is, the 
NI transition temperature is found to decrease (see Table I). Therefore, when comparing 
b for different values of is, it is appropriate to rescale the temperature T* and compare 
measurements with comparable values of the reduced temperature r (resulting in similar 
values of s&). 

To facilitate the analysis, we decided to remove the wall at z = allowing us to deal 
with intrinsic anchoring alone. When v % 0.3 (but nonzero) the intrinsic easy axis remains 
homeotropic [11], and the hybrid cell-like geometry studied in the previous section is main- 
tained. As in the case of pure external anchoring, we find that the intrinsic extrapolation 
length b shows no temperature dependence (within estimated error; see Fig. 4). However, b 
does depend on the anisotropy of the nematic-nematic interparticle potential: for v = 0.05 
and v = 0.1 one finds b = 16.3a(l ± 13%) and b = 4.6a(l ± 10%), respectively (see Table I 
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FIG. 5: Same as Fig. 4, but for different values of the bulk interaction anisotropy: v = 0, v = 0.05, 
and v = 0.1 (diamonds, squares and triangles, respectively). The corresponding temperatures are 
T* = 1.34, T* = 1.195, and T* = 1.065, respectively, ensuring that the bulk value of the order 
parameter and the reduced temperature be the same in all cases (r ~ 0.12). There is no intrinsic 
anchoring for u = 0, while for v ^ an increase of v results in a decrease of the extrapolation 
length. 

and Fig. 5). This trend can be attributed to an increase of the anchoring energy W, as well 
as to a simultaneous decrease of the elastic constant K upon increasing v. Note that in 
the analysis of Ref. [11] b = 8a(l ± 13%) was found for v — 0.1 at zero temperature. The 
disagreement with the current largely temperature-independent estimate can be attributed 
to considering only nearest neighbors in the present analysis, which — as already stated - 
results in an underestimation of b. 

Note, moreover, that for v > 0.1 extrapolation lengths are in the microscopic range (of 
the order of a few — up to 5 — lattice spacings a) . On the other hand, experimental values 
of the extrapolation length are typically of the order of 100 nm or greater [9]. We can obtain 
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quantitative agreement between our results and experiments by a significant decrease of the 
v parameter, as also suggested in Ref. [11]. A small value of v in Eq. (1) promotes parallel 
alignment, as is favored, e.g., by steric repulsions in a system of hard rods. A decrease in v 
might therefore be regarded as an effective inclusion of steric repulsions which are absent in 
our lattice model. 

V. CONCLUSIONS 

In this paper we have studied the anchoring of a nematic liquid crystal to a solid substrate 
and to a free interface, using numerical simulations of a simple hexagonal lattice model of 
two-dimensional spins interacting through a spatially anisotropic potential. We focused on 
the roles of substrate roughness and of the spatial anisotropy of the interparticle potential. 
An elastic deformation was imposed on the simulation cell, allowing us to extract the surface 
extrapolation length b from the director profile. Setting the anisotropy of the potential 
to zero allowed us to probe the strength of the external anchoring (the anchoring arising 
from the direct interaction of the nematic with the surface). For temperatures not too 
close to the nematic-isotropic transition the extrapolation length was found to be largely 
temperature independent, and to grow with increasing surface roughness (as characterized 
by the distribution of the local preferred axis on the surface). This increase is physically 
reasonable, reflecting the decrease in the overall anchoring of the nematic as the surface 
becomes more disordered. Qualitatively similar behavior was found when we simulated the 
original Lebwohl-Lasher model where the spins are three-dimensional. 

With the full anisotropic interaction potential present, intrinsic anchoring arises due to 
the incomplete spin-spin interactions at the sample surface. For a free nematic interface 
(and external anchoring absent), we extracted the intrinsic surface extrapolation length. 
Like the extrapolation length in the presence of only external anchoring, this length is 
temperature-independent. We found that the strength of the intrinsic anchoring grows 
together with the elastic softness of the nematic as the interaction potential anisotropy is 
increased, leading to smaller values of the extrapolation length. Thus, obtaining agreement 
with the experimentally measured values of b, on the order of 100 nm, requires a relatively 
small value of the interaction anisotropy parameter v approaching 10~ 2 . 
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